rm(list = ls())

library(lfe)
library(texreg)
library(here)
library(data.table)
library(qs)
library(tidyverse)

library(splm)
library(spdep)
library(plm)

# Load and preparedata ---------------------------------------------------------------
#Voting and immigration by municipality
cb <- qread(here("data", "cb.qs")); setDT(cb)

#Cross walk data to match new with early municipalities
crosswalk <- read_csv(here("data", "cw_munnr12_bfs19.csv"))

#Coordinate data for early municipalities
coords <- haven::read_dta(here("data", "cleaning", "auxiliary_files", "a6_munnr2012_coordinates.dta"))

#Add bfs19 identifier
coords <- coords %>% left_join(crosswalk)

# Calculate new centroids at bfs19 level
centroid_df <- coords %>%
  group_by(bfs19) %>%
  summarise(
    longitude = mean(longitude),
    latitude = mean(latitude))

#Add centroid data to cb
cb <- cb %>% left_join(centroid_df)

#Filter to election years
cb <- cb %>% filter(complete.cases(frPrct))


#Filter to include only border region and travel distance less than 30 minutes
cb_filtered <- cb[BR == 1 & travelMUNmin <= 30,]


# Prepare spatial weights -------------------------------------------------
#Find unique municipalities to eliminate duplicate entries for the same location
unique_municipalities <- cb_filtered %>% dplyr::select(bfs19, longitude, latitude) %>% distinct()

#Extract longitude and latitude
coordinates <- unique_municipalities[, .(longitude, latitude)] 

#find the k-nearest neighbors for each municipality based on their geographic coordinates
knn_weights <- knearneigh(coordinates, k = 4) 

#If error (because of different SPLM version), uncomment the following lines
#coordinates <- cbind(as.numeric(coordinates$longitude), as.numeric(coordinates$latitude))
#knn_weights <- knearneigh(coordinates, k = 4) 

#Convert to a neighobor-list
nb <- knn2nb(knn_weights) 

#Transform into spatial weights
listw <- nb2listw(nb, style = "W", zero.policy = TRUE) 


# Regressions accounting for spatial dependence -------------------------------------------------

#Spatial Panel Fixed Effects Error Model, estimated via Maximum Likelihood. Addresses potential spatial autocorrelation in the error terms.
fmla <- frPrct ~ transBorder15 + freeBorder15


#spatially lagged dependent variable model
mod_spml1 <- spml(fmla, data = cb_filtered, listw = listw, model="within", effect="twoways",lag=TRUE,
                  spatial.error="none") 
coefficients1 <- coef(mod_spml1)


#Sptatial autocorrelation in error term
mod_spml2 <- spml(fmla, data = cb_filtered, listw = listw, model="within", effect="twoways",lag=F,
                 spatial.error="b") 

# Extract coefficients, standard errors, and p-values from both models
coefficients1 <- coef(mod_spml1)
standard_errors1 <- sqrt(diag(vcov(mod_spml1)))
p_values1 <- 2 * (1 - pnorm(abs(coefficients1 / standard_errors1)))

coefficients2 <- coef(mod_spml2)
standard_errors2 <- sqrt(diag(vcov(mod_spml2)))
p_values2 <- 2 * (1 - pnorm(abs(coefficients2 / standard_errors2)))

# Create texreg objects for each model
texreg_obj1 <- createTexreg(coef.names = names(coefficients1), coef = coefficients1, se = standard_errors1, pvalues = p_values1)
texreg_obj2 <- createTexreg(coef.names = names(coefficients2), coef = coefficients2, se = standard_errors2, pvalues = p_values2)

# Print the regression tables
screenreg(list(texreg_obj1, texreg_obj2))


source("functions.R")
Texreg(list(texreg_obj1, texreg_obj2), custom.header = list("Percent Support for Anti-Immigrant Parties" = 1:2), 
       custom.gof.rows = list("Model" = c("Spatial Correlation in Error Term", "Spatially Lagged DV")),
       custom.coef.map = list("rho" = "\\rho", "lambda" = "\\lambda", "transBorder15" = "0--15 Minutes X Transition", "freeBorder15" = "0--15 Minutes X Free"),
       caption = "Effects on voting, accounting for spatial dependence.  Outcome: Percent support for anti-immigrant parties in federal parliamentary elections. Sample: Border region, travel minutes $<=$ 30, 1991--2019 elections.",
       label = "tab:spatial_dependence",
       file = here("tables", "spatial_dependence.tex")
       )






